Lifetime of double occupancies in the Fermi-Hubbard model 
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We investigate the decay of artificially created double occupancies in a repulsive Fermi-Hubbard 
system in the strongly interacting limit using diagrammatic many-body theory and experiments with 
ultracold Fermions on optical lattices. The lifetime of the doublons is found to scale exponentially 
with the ratio of the on-site repulsion to the bandwidth. We show that the dominant decay process 
in presence of background holes is the excitation of a large number of particle hole pairs to absorb 
the energy of the doublon. We also show that the strongly interacting nature of the background 
state is crucial in obtaining the correct estimate of the doublon lifetime in these systems. The 
theoretical estimates and the experimental data are in fair quantitative agreement. 

PACS numbers: 03.75.Ss, 05.30.Fk, 34.50.-s, 71.10.Fd 



The non-equilibrium dynamics of a strongly interact- 
ing quantum many-body system is one of the most com- 
plex problems of modern physics. It encompasses var- 
ious fields from the cosmology of the early universe [T] 
or non-equilibrium jet production in high energy heavy 
ion collisions [2] to pump-probe experiments and opera- 
tion of solid state devices under strong drive [3] in con- 
densed matter physics. There are many open questions 
concerning non-equilibrium processes from both a theo- 
retical and an experimental perspective, especially in the 
realm of condensed matter physics. 

The theoretical understanding of interacting quantum 
many body systems in thermal equilibrium is on a much 
stronger footing, although strongly interacting systems 
like high temperature superconductors are not yet com- 
pletely understood. This understanding is based on 
paradigms such as the quasiparticle excitations in the 
Fermi liquid model and ground states with broken sym- 
metry described in terms of order-parameters and their 
fluctuations. The crucial point in all these paradigms 
is the hierarchy of energy scales of the quantum states. 
By working with a restricted set of states, organized ac- 
cording to their energy, it is possible to obtain a simpli- 
fied model of the system. These low energy descriptions 
can capture the response of the system under small per- 
turbations from equilibrium. However, in systems far 
from equilibrium, there is no organizing principle as the 
dynamics couples disparate states with widely different 
energies and linear response theory breaks down. This 
makes it hard to construct generic paradigms and one 
needs to solve the full microscopic Hamiltonian dynam- 
ics of an interacting quantum many-body system. 

Some progress has been made for ID systems, where it 
is often possible to obtain exact solutions for the eigen- 
states of the Hamiltonian. The absence of thermaliza- 
tion in ID Bose systems has been predicted [H |5] and 
observed [5] in cold atomic gases. However, these studies 



are hard to generalize to higher dimensions. 
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FIG. 1: Stability of highly excited states in the single-band 
Hubbard model. Doubly occupied lattice sites are protected 
against decay by the on-site interaction energy U . The aver- 
age kinetic energy of a single particle in a periodic potential 
is half the bandwidth 6 J. Thus the relaxation of excitations 
requires several scattering partners to maintain energy con- 
servation. 

In this context, it is useful to seek answers to con- 
crete and focused questions involving non-equilibrium 
dynamics of specific strongly interacting systems. They 
have practical importance and help us gain better un- 
derstanding of classes of non-equilibrium processes. Re- 
cent advances in controlling ultracold atomic gases with 
and without optical lattices have led to their emergence 
as perfect systems to study such phenomena. These 
systems, which can simulate strongly interacting model 
Hamiltonians, are essentially decoupled from external 
heat baths and hence the intrinsic non-equilibrium dy- 
namics of the system can be studied easily. Compared to 
condensed matter systems, the low density in these sys- 
tems results in long timescales for dynamics. As a result 
the system can be followed in real time without the use 
of ultrafast probes. Further, it is relatively easy to cre- 
ate and characterize an initial state far from the ground 
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state, which is crucial since the dynamics depends heavily 
on the initial state. 

In fact, questions of non-equilibrium dynamics and 
thermalization timescales are particularly important for 
these artificially engineered strongly correlated systems. 
Their key feature is the precise tunability of the Hamilto- 
nian parameters which has made these systems ideal for 
the simulation of strongly interacting many-body Hamil- 
tonians relevant to condensed matter systems. However, 
an implicit assumption in this comparison is that the 
system is in thermal equilibrium at low temperatures. In 
this context it is important to estimate the thermaliza- 
tion timescales as these systems are always characterized 
by a finite sample lifetime. Besides, several proposed 
methods to prepare the system in novel phases explicitly 
depend on adiabatic tuning of Hamiltonian parameters, 
which place stronger constraints on the possible sweep 
rates than mere demand of thermalization. 

An important class of non-equilibrium problems is the 
decay of a high energy excitation into low energy excita- 
tions. This problem occurs in diverse contexts like multi- 
phonon decay of excitons in semiconductors [7 , pump 
and probe experiments [3^ and dynamics of nuclear res- 
onances [5^. In this paper, we study this problem in the 
non-equilibrium dynamics of artificially created double 
occupancies in the Fermi Hubbard Model in the strongly 
interacting regime. Specifically, we will look at the mech- 
anism of doublon decay in this system and the relation 
of the doublon lifetime to the repulsive interaction. We 
study this dynamics both experimentally using ultracold 
Fermions on an optical lattice [9] and theoretically using 
a projected Fermion model and diagrammatic resumma- 
tions. 

The doublon lifetime has practical implications for 
the sweep rates of Hamiltonian parameters in cold atom 
systems in the following way: The usual access to the 
strongly interacting regime is to start with a weakly in- 
teracting system and increase the ratio of interaction U 
to the hopping energy J. As this ratio increases, the 
density of doublons in the system in equilibrium should 
decrease. Thus the doublon lifetime provides the dom- 
inant equilibration timescale for the system. We note 
here that this problem has structural similarities with the 
decay of a deeply bound excitonic state through multi- 
phonon processes in semi-conductors [7], but as we shall 
see, the strong Hubbard repulsion modifies the situation 
in an essential way. 

Our main results are (i) The decay of a doublon is a 
slow process as the doublon needs to distribute a large en- 
ergy (~ U) to other excitations in the system which have 
a much smaller energy scale, (ii) The primary mode of de- 
cay of the doublon involves creation of particle-hole pairs 
in the background system, (iii) The decay rate scales as 
T ^ C J exp(— a U/J) and the decay becomes slower with 
increasing interaction. We obtain C and a from experi- 
ments and from theoretical calculation, (iv) We find that 
the interactions between pairs of single Fermions, which 
in our model are induced by projection, are important 



and quantitatively affect the timescale of the decay. Thus 
the strongly correlated character of the system changes 
the dynamics in an essential way. 

The paper is organized as follows: In Section I we dis- 
cuss the various possible decay mechanisms of the dou- 
blon in these systems and give a scaling argument for the 
decay rate in each case. In Section II we describe the ex- 
periments and its results. In Section III we discuss the 
most relevant decay mechanism in our experiments and 
develop the theoretical model for doublon decay. In Sec- 
tion IV we outline the diagrammatic method to compute 
the doublon lifetime. In Section V we discuss the theoret- 
ical results and its comparison with the experiments. We 
conclude in section VI by discussing the importance of 
these results and future directions. The technical details 
of the theory are described in relevant appendices. 

I. DECAY MECHANISMS FOR A DOUBLON 

The single-band Hubbard model describing the 
Fermions on an optical lattice is given by |15| 

H = -JY^ cl^Cj^ + UY^ n.^m^. (I) 

At large C//J, this model has three main energy scales. 
There is the energy of double occupancies, given by the 
Hubbard repulsion [/, the kinetic energy of the Fermions 
given by the tunneling J and the superexchange scale 
Jea; — AJ'^ /U , which governs the spin dynamics in the 
system. At large C//J, these scales are well separated 
from each other, U ^ J ^ J^x- As we show below, the 
separation of the energy scale U from the other energy 
scales J and Jex leads to a slow decay of doublons in the 
system. 

In order to decay the doublon has to give up its energy 
^ U to other excitations in the system. Let the typical 
energy of a possible excitation be eo where eg can be 
either ~ J or ~ J^x depending on the background state 
in which the doublon is propagating. We assume that 
eo ^ U, so that a large number n ^ U /eo of excitations 
must be created to satisfy energy constraints. The matrix 
element for this process can be calculated by an n*'' order 
perturbation theory and is given by 



The decay rate is ^ in units of J . Using Stirling's 
formula, and the fact that neo = t/, we find that for large 
n the decay rate scales as 

r - J{ J/U)^ - CJexp(-a U/eo log{U/J)) (3) 

where C and a are constants which we will extract from 
detailed calculations and experimental data. 

In order to discuss the specific decay mechanisms of a 
doublon, we need to specify the state of the background 
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system in which the doublon is propagating. If the sys- 
tem is a homogeneous Mott insulator at half-fiUing, the 
only possible candidate for transfer of energies are spin 
excitations with bandwidth eo ^ J ex- This leads to the 
decay rate scaling as F ^ Jexp(— a U"^ j ,P log(f7/ J)) and 
is an extremely slow process. However, if the system is 
compressible, the dominant energy transfers are to ki- 
netic energy of the Fermions through creation of particle- 
hole pairs with typical energy eg ~ J. This leads to the 
decay rate scaling as F ~ Jexp(— a ?7/Jlog(C//J)). This 
is a much faster decay process and will dominate over de- 
cay through spin excitations. We note that compressible 
states with holes can exist (i) at the edges of systems with 
confining traps or (ii) in the bulk of the system as a re- 
sult of a large density of doublons created by modulation 
spectroscopy. In a trapped system, there is another pos- 
sibility of giving up the energy to the potential energy of 
the Fermions at the edges. This however involves trans- 
fer of particles from the center to the edges of the trap 
and is usually a much slower process for typical shallow 
traps used in cold atom experiments. 

As we will see in the next section our experimental 
system has a lot of holes and we can eliminate many 
of the possible decay mechanisms for our experimental 
configuration. Therefore, in this Paper we shall focus on 
the dominant doublon decay channel involving excitation 
of particle hole pairs. 



II. EXPERIMENTS 

This section describes the experimental steps towards 
the observation of doublon relaxation: Initially, a sam- 
ple of repulsively interacting, ultracold fermions is pro- 
duced and loaded into an optical lattice. Starting from 
this equilibrium state, we create additional double occu- 
pancies via lattice modulation. Immediately after this 
excitation we measure the time evolution of the dou- 
ble occupancy. We remove the influence of inelastic loss 
processes by comparing to a reference measurement and 
extract the elastic doublon lifetime using a simple rate 
equations model. Finally, this elastic lifetime is normal- 
ized with the tunneling time Jjh and found to depend 
exponentially on J7/6J. 

The experimental sequence used to produce quantum 
degenerate Fermi gases has been described in detail in 
previous work In brief, we prepare (50 ± 5) x 10"^ 

''^K atoms at temperatures below 15% of the Fermi tem- 
perature Tp in a- balanced mixture of two magnetic sub- 
levels of the F = 9/2 hyperfine manifold. The confine- 
ment is given by a crossed beam dipole trap with trap- 
ping frequencies i-^x,y,z = 27r x (35, 23, 120) Hz. To ac- 
cess a wide range of repulsive interactions we make use 
of two magnetic Feshbach resonances. With mp = 
(—9/2, —1/2) spin mixture, we realize scattering lengths 
of 98 ao and 131 ao, where ag is the Bohr radius [12]. The 
(—9/2, —5/2) spin mixture allows us to reach the strongly 
repulsive regime with scattering lengths of 374 oq, 571 oq 



and 672 ao After adjusting the scattering length 

to the desired value, we add a three-dimensional optical 
lattice of simple cubic symmetry. The lattice depth is 
increased in 200 ms to final values between Q.hEpi and 
12.5i?ij in units of the recoil energy Ej^ = h?/2m\^. 
Here h is Planck's constant, m the atomic mass and 
A = 1064 nm the wavelength of the lattice beams. The 
lattice beams have Gaussian profiles with 1/e^ radii of 
Wx,y,z = (160, 180, 160) iim. at the position of the atoms. 
For a given scattering length and lattice depth, J and U 
are inferred from Wannier functions [TT1[TS]. Their statis- 
tical and systematic errors are dominated by the lattice 
calibration and the accuracies in width and position of 
the two Feshbach resonances [121 IE] ■ 

Depending on U and J the final states of the system 
range from metallic to Mott insulating phases, but always 
with a double occupancy below 15%. This equilibrium 
system is now excited by a sinusoidal modulation of the 
lattice depth with a frequency close to U/h. This causes 
an increase of the double occupancy between 5 and 20% 
as compared to the initial state. The modulation am- 
plitude is 10% on all three axes, while the modulation 
duration was chosen such that the amount of doubly oc- 
cupied lattice sites saturates [HI [T6HI9] . The system is 
now in a highly excited non-equilibrium state with double 
occupancies between 15 and 35%. 

After free evolution at the initial lattice depth and in- 
teraction strength for a variable hold time up to 4 s we 
probe the remaining double occupancy of the system. 
This is accomplished by a sudden increase of the lattice 
depth to 30£'r, which prevents further tunnehng. We 
then measure the amount of atoms residing on singly 
(doubly) occupied sites TVs (A^d) by encoding the double 
occupancy into a previously unpopulated spin state us- 
ing radio frequency spectroscopy [TT. Combining Stern- 
Gerlach separation and absorption imaging we obtain 
the single occupancy Us = Ns/Ntot, double occupancy 
rid — Ad/Atot and total atom number A^tot = A^s + Ad- 

The time evolution of the double and single occupancy 
and of the total atom number is shown for two differ- 
ent parameter sets in the upper row of Fig. [2] In both 
cases, the double occupancy decays exponentially within 
the observation time, and the single occupancy rises ac- 
cordingly. The time evolution of the total atom num- 
ber, however, exhibits a remarkable difference between 
the (mp = -9/2, -7/2) and the (mp = -9/2, -5/2) spin 
mixture: Whilst the atom number of the (-9/2,-7/2) 
sample remains rather constant, the (—9/2, —5/2) sample 
suffers from an atom number reduction of 50 % within 2 s. 
This behavior can be observed for all parameter sets and 
is a consequence of the shorter lifetime of a (—9/2, —5/2) 
spin mixture. 

The only relevant process described by the Fermi- 
Hubbard model is the decay of a doublon into two single 
particles which remain within the system. The time as- 
sociated with this process will be called doublon lifetime. 
In an experiment, inelastic processes may occur, result- 
ing in atoms exiting the system. For a valid comparison 
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FIG. 2: Time evolution of double occupancy, single occupancy and total atom number for different ratios U /6 J. In the upper 
row, the system was previously excited via lattice modulation. The bottom row shows the reference measurement for the 
determination of the residual dynamics. The round data points were recorded using a mp = (—9/2, —7/2) spin mixture with 
U/h — 1.4kIIz and J/h — 70 Hz, whereas the triangular data points show a (-9/2,-5/2) mixture with U/h = 3.2kHz and 
J/h — 100 Hz. The solid lines are simultaneous fits of the integrated population equations of Eq.|4] The total atom numbers are 
scaled to the initial values. Single occupancy and double occupancy are the fraction of atoms residing on sites of the respective 
type. Due to different detection efficiencies for hyperfine states the sum of double and single occupancy can be higher than 
one. Error bars denote the statistical error of at least four identical measurements. 



with theory it is therefore crucial that these processes 
do not interfere with the determination of the doublon 
lifetime. In the following, we show how we eliminate the 
influence of inelastic loss processes on the observation of 
the doublon decay. 

For every dataset on doublon decay after lattice modu- 
lation, we record a corresponding reference dataset with- 
out lattice modulation, but with the same system pa- 
rameters. Two of these reference datasets are presented 
in the bottom row of Fig. [2] They show the dynamics 
of double occupancies and atom number governed by in- 
elastic processes, which are not taken into account by the 
Fermi-Hubbard model. 

Combining these two measurements, we can unambigu- 
ously extract the doublon lifetime by simultaneously fit- 
ting a system of coupled rate equations. They describe 
the population dynamics in the optical lattice, consider- 
ing three general processes: 
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The total number of atoms on doubly occupied sites N^i 
is written as the sum of the equilibrium population N^.q 
and the additional amount of double occupancy AN^ cre- 
ated by the lattice modulation. The three time constants 
correspond to three independent local decay processes 
differing in the type of site they affect: td describes the 
population flow from doubly occupied to singly occupied 
lattice sites visible as a decay of double occupancy within 
0.01 — Is that is accompanied by a rise of the single oc- 
cupancy. We identify this time with the lifetime of dou- 
blons. The other two times denote loss time constants, 
which lead to a reduction of the total atom number: rioss 
corresponds to losses affecting both site types in the same 
manner, which is only observed in the total atom num- 
ber. Additional inelastic losses on doubly occupied sites 
are summarized by Tin, visible as a simultaneous decay of 
both the total atom number and double occupancy. This 
model does not account for changes of the decay times 
during the decay or for higher order terms in the rate 
equations. 

Since the modulation has no influence on the evolu- 
tion of the total atom number, this procedure removes 
the influence of Tin and rioss • A reliable determination of 
the doublon lifetime td is thus possible if it differs sig- 
nificantly from the loss times. The model and the obser- 
vation are found to agree very well within experimental 
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uncertainties, as can be seen in Fig. [2j 

We measure this doublon lifetime for various tunnel- 
ing and interaction strengths, covering a parameter range 
where J and U each differ by at least a factor of four. 
The determined lifetimes vary over two orders of mag- 
nitude, as shown in Fig. [3j Furthermore, the lifetime 
clearly does not depend on the tunneling energy or the 
interaction energy alone. 
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FIG. 3: Doublon lifetime as a function of U and J. The 
round data points show the fit results to datasets as shown 
in Fig .[2] obtained with a (—9/2,-7/2) spin mixture while 
the triangular points correspond to the (—9/2, —5/2) mixture. 
Error bars denote the confidence intervals of the lifetime fits 
and the statistical errors in U and J. 

Since the tunneling time h/ J is the dominant timescale 
of dynamics in an optical lattice, it appears natural to 
express the doublon lifetime in units of h/ J . After this 
rescaling, we found that, to a good extent, the doublon 
lifetime only depends on the ratio U /&J. 

Fig.|4]shows the doublon lifetime in units of the tunnel- 
ing time versus U /6 J on a logarithmic scale. Remarkably, 
over the entire parameter range the data collapses in a 
corridor and can be described by an exponential function 
of the form: 



^ = C exp [a^^ 



(5) 



The scaling exponent a is found to be a = 0.82 ± 0.08 
with C — 1.6 ± 0.9. This is in reasonable quantitative 
agreement with the following calculation of the doublon 
lifetime. 

The slight offset between the two spin mixtures in 
Fig. |4] could be due to the fact that the absolute values 
for U and J differ significantly between the (—9/2, —5/2) 



and the (-9/2,-7/2) mixture [20]. Whilst the ratio be- 
tween interaction energy and kinetic energy J7/6 J, which 
dominates the dynamics, lies in the same range, the ab- 
solute values also matter in an inhomogeneous system. 
For the (-9/2,-7/2) mixture the higher ratio of chem- 
ical potential to on-site interaction is expected to lead 
to a higher filling in the trap centre and consequently 
to a higher equilibrium double occupancy N^fi than for 
the (—9/2, —5/2) mixture. It is conceivable that this dif- 
ference modifies the dynamics of doublon creation and 
doublon relaxation. 

In additional measurements we examined the depen- 
dence of the doublon lifetime on the initial double occu- 
pancy AA^d and on the total atom number N . In the 
former case, we reduced the lattice modulation ampli- 
tude from 10% to 5%, resulting in AA^^ — 9% instead 
of AA^rf = 18%, while keeping all other parameters con- 
stant with U/6J = 4.5. The measured lifetimes agree 
within the error bars, they are td_5% = (77 ± 25) x h/J 
and td 10% = (58 ± 10) x h/J, respectively. In the lat- 
ter case, we prepared two otherwise identical samples 
at U/6J = 3.4 with A^ = (49 ± 7) x 10^ atoms and 
with A^ = (26 ± 4) X 10'^ atoms, respectively, yielding 
TD,490oo = (11 ± 2) X h/J and rD,26000 = (19 ± 2) x h/J. 

This shows that, although there is a dependence on the 
total density and on the doublon density, these effects 
are small compared to the dominant scaling with U /6 J. 
Their systematic study is beyond the scope of this work. 



III. THEORETICAL MODEL OF DOUBLON 
DECAY 



We consider the decay of an isolated doublon moving in 
the homogeneous background of a compressible state of 
single Fermions. Before constructing a model for doublon 
decay, we focus on the dominant mechanism of decay. In 
the experiments, lattice modulation created 15 — 35% 
double occupancies. Assuming an initial half-filled sys- 
tem, half the amount of holes were also created in the 
system. At these hole densities, the kinetic energy as- 
sisted decay scaling as ~ exp(— 17 /J) is much faster than 
the spin fluctuation or doublon-doublon collision assisted 
decay which scale as ^ exp(— J7^/J^) [10_. Further, the 
population of higher bands can be excluded, since U is 
always smaller than half the band gap. We also note that 
as the diff'erence between U and the chemical potential is 
always positive, confinement assisted decay of doublons 
near the edge of the cloud is unlikely, as the accessible 
confinement energy is not very large, and the tunneling 
rate is very small. Finally, a homogeneous compressible 
background is justified since most of the doublons are 
created in the central region of the trap, where the filling 
is highest, and decay at most within a few sites of where 
they are produced The estimated travel distance for a 
random walk during the decay process is not more than 
y^T^J/h < 10 sites, which is less than the cloud radius. 

In our experiments, the doublons and holes are created 
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at finite density by driving the system with optical lattice 
modulations. The relaxation of the system to equilibrium 
involves two very different time scales. The first timescale 
is associated with the relaxation of holes and doublons 
to a state of quasi-equilibrium without the decay of dou- 
blons. The second timescale, which is the focus of this 
paper, is associated with the decay of doublons into sin- 
gles. We expect that the second timescale is much slower 
than the first. Moreover, we expect that non-linear ef- 
fects of doublon decay as doublon-doublon scattering can 
be neglected since their kinetic energy ~ ,P /U is small. 
Thus in this paper we consider the problem of the de- 
cay of a single doublon in the background of equilibrated 
Fermions. 

To construct our model Hamiltonian, we explicitly 
treat the doublon as a separate entity from the back- 
ground Fermions. This approximation is justified in the 
strongly interacting limit due to the separation of dou- 
blon and background Fermion time-scales. We split the 
complete Hamiltonian of the system into three parts 



H = Hf + Hd + H 



(6) 



where Hf describes the background Fermion 
subsystem — which we model as the projected Fermi 
sea, Hd describes the on-site interaction of the pair of 
Fermions that make up a doublon, and Hfd describes 
the Fcrmion-doublon interaction. The details of how 
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FIG. 4: Semilogarithmic plot of doublon lifetime td vs. 
U/6J. The lifetime is extracted from datasets as shown in 
Fig. PI Solid and hollow circles denote the (—9/2,-5/2) and 
(—9/2, —7/2) spin mixture respectively, while the dashed line 
shows the theoretical result at half filling. The solid line is a fit 
of Eq. [5] to the experimental data, yielding a — 0.82 ± 0.08, 
whereas for the theory curve the asymptotic slope at large 
U/6J is qt = 0.80. The shaded corridor was obtained by 
varying the filling factor in the calculation by 0.3. This has 
only a weak effect on the slope. The inset shows the param- 
eters used to realize the different values of U /6 J. Error bars 
denote the confidence intervals of the lifetime fits and the sta- 
tistical errors in U /6 J. The systematic errors in U /6 J and 
Td = h/J are estimated to be 30% and 25%, respectively. 



to separate the Fermi-Hubbard Hamiltonian into the 
above three parts via projection operators are discussed 
in Appendix |^ The projection operators induce inter- 
actions in the Fermion subsystem as well as between 
the Fermions and the doublons. The Fermion doublon 
interactions are responsible for the doublon decay, and 
the Fermion-Fermion interactions modify the lifetime 
substantially. 

As mentioned above, we expect hole density in these 
systems to be ~ 15%. At such high hole densities the 
projected Fermi sea is a good approximation for the back- 
ground state. Further the temperature of the system is 
high enough (T J) to prevent formation of more or- 
dered states like superfluids. 

Except for the single doublon that is undergoing de- 
cay, the large energy cost of double occupancies is taken 
into account by projecting out configurations with dou- 
ble occupancies from a simple Fermi sea. In the pro- 
jected subspace, the Fermions can only hop in the pres- 
ence of empty sites (holes) and are governed by the ef- 
fective Hamiltonian 

Hf = -J ^{'^-n^cr)c\„c.J„{l~nj^)^^l'^c\„c^^ (7) 

where creates a Fermion with spin cr, n^o- is the 
corresponding number operator, and /i is the chemi- 
cal potential. Expanding out the Hamiltonian one gets 
Hf 



H^ + Hp, with 



H° - E 4s^-mE4c.., (8) 

{ij),a- i,cr 

Hp = JiY^ rii^cl^Cj^ + cl^c^^Uj^, (9) 

where we have replaced J by Ji in the second term. Ji 
will be treated as a perturbation parameter to organize 
the calculation but we will put Ji = J at the end of the 
calculation. Hp, coming from the projection operators 
can thus be interpreted as a Fermion-Fermion scattering 
term which leads to the creation of particle-hole pairs. 
We thus see that projection induces interaction between 
the Fermions. 

We note that the scattering is always between Fermions 
of opposite spins. Since we will be interested in calcu- 
lating Feynman diagrams, we note that the interaction 
vertex for the Fermion Fermion scattering can be written 
as Ji (7k + 7k-q) , where k and k — q are momenta of the 
incoming and outgoing Fermion with the same spin and 

7k — 2 X]t=i i'^ ^ dimensions. This is depicted in 

first row of Table H] 

Throughout our treatment, we leave out terms such as 
— J ni-ffcl^Cj^rij-^ in Eq. ^ involving six or more fermion 
creation/annihilation operators. Intuitively such terms 
are rare because they involve collisions of multiple parti- 
cles. 

We now consider the decay of a single doublon in this 
background state. The doublon (d) and Fermion-doublon 
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Interaction 



Fermion-Fermion 
Scattering 



Doublon-Fermion 
Scattering 



Doublon Decay 



Diagram 



iH.a Pi + ^-1 - k-2., a 



p-q + k 




p - k. I 



Vertex 
Momentum 
Avg. Vertex 



J (7ki + 7k2 



J(7p-q +7p +7q) 



J(7k + 7p-k) 



TABLE I: Interaction vertices for different processes in the 
model for doublon decay. The single lines are Fermion prop- 
agators while the double lines are doublon propagators. The 
top entry in the right-most column is the corresponding ver- 
tex function, while the bottom entry is the "momentum av- 
eraged" vertex function that we use in the our resummation 
technique. The first row corresponds to Hp (Eq. [9| while the 
next two rows correspond to the two terms that make up Hfd 
(Eq. 111. Here yk ~ 2{cos{kx) + cos(fcy) -I- cos(A;z)) and z is 
the coordination number, z — 6 for 3D cubic lattice. 



i fd) Hamiltonians can be written as 

i 

Hf, = jj2{4d^+d]d^ + d]d 



ft r . 



(10) 

(11) 



where H.c. stands for the Hermitian conjugate of the pre- 
ceding term. The doublon interacts with the Fermions in 
two different ways: (i) it can scatter off a Fermion leading 
to the hopping of doublons with back-flow of Fermions; 
(ii) it can decay by creating a singlet particle-particle 
pair. 

The interaction vertices of the doublon with the 
Fermions are given in second and third rows of Ta- 
ble |l] The vertex for scattering off particle-hole pairs 
is J(7p_q -l- 7q -l- 7k), where p is the momentum of the 
incoming doublon, q is the momentum of the outgoing 
Fermion and k the momentum of the outgoing hole. The 
corresponding vertex for decay through singlet creation 
is given by J(7k + 7p-k) where k and p — k are the 
momenta of the Fermions created. 

We assume that we are looking at the decay of a single 
doublon i.e. while the doublon is affected by the pres- 
ence of the background Fermions, the Fermions are un- 




FIG. 5: A Typical doublon self-energy diagram. The double 
lines are bare doublon propagators while the single lines are 
bare Fermion propagators. The dashed line cuts the diagram 
in half and shows the final products of the process represented 
by this diagram 



affected by the presence of the doublon. The motivation 
for this assumption is that the experimentally observed 
decay rate depends only weakly on the doublon density. 



IV. DIAGRAMMATIC COMPUTATION OF 
DOUBLON LIFETIME 

Our strategy for finding the lifetime of a doublon is to 
calculate the doublon Green function 



(12) 



where is the self-energy arising from interaction with 
Fermions. The imaginary part of the self-energy at 
oj = U then gives the decay rate F and its inverse is 
the required lifetime r. Since we are interested in the 
high frequency response, the momentum dependence of 
the self energy should be negligible in this limit. 

We perform the calculation at T = 0, where the rela- 
tion between imaginary part of the self-energy and decay 
rate is exact. At finite temperatures Im S(a;) has an ex- 
tra contribution due to scattering on particle-hole pairs 
created by thermal fluctuations. Thus, we must com- 
pute the scattering rate separately, and subtract it from 
Im S(w) to obtain the decay rate. However, since we are 
looking at frequencies ~ U , ignoring thermal fluctuations 
is justified for T , which is the regime of interest. 

Physically, there are two important processes for the 
doublon decay. A doublon can lose its energy either by 
creating a large number of particle-hole pairs, each with 
an energy ~ J, or by creating a few high energy particle- 
hole pairs, each of which is unstable and creates a shower 
of particle hole pairs of low energies. The first process 
is a high order diagram in the doublon self-energy while 
the second process comes from high order diagrams in the 
Fermion self-energy. We find that combinations of both 
processes give important contributions to the doublon 
decay rate. 

The typical doublon self-energy diagram (Fig. [5| de- 
picts a process of creation of a number of particle and 
hole excitations. Since we are interested in the imagi- 
nary part of the self energy at uj = U, the Fermion lines 
crossing the dashed line which cuts the diagram in half 
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FIG. 6: Self-consistency equation for the doublon propagator 
(top) and some typical diagrams that make up the full prop- 
agator (bottom). Thin double-lines indicate the bare dou- 
blon propagator, the double-lines with a squiggle the full (re- 
summed) doublon propagator, and thick single-lines the full 
(resummed) Fermion propagator. 



should be on-shell and their energies must add up to U. 
The leading order contributions to the decay rate thus 
come from the diagrams which maximize the number of 
Fcrmions that cross the dashed line while minimizing the 
number of interaction vertices. 

Our approach for obtaining the doublon self energy 
consists of (1) obtaining the projected Fermi sea Green 
function, and (2) using it to obtain the doublon self- 
energy. We make the dilute doublon approximation, and 
assume that the Fermion Green function is independent 
of the doublon Green function. We proceed by formulat- 
ing a diagrammatic resummation technique for the dou- 
blon self-energy in the following subsection. In doing so, 
we relate the doublon self-energy to the Fermion Green 
function, which we calculate in the next two subsections. 



A. Doublon Self-Energy 

For large U/ J, doublon decays into a large number 
of particle-hole pairs, and therefore one needs to com- 
pute high order diagrams to obtain the doublon self- 
energy (for creation of n pairs, one needs to compute 

2n! diagrams). It is then much preferable to resum 
a class of diagrams, rather than evaluate an exponen- 
tially increasing number of them. We use a self-consistent 
non-crossing approximation to achieve this resummation. 
The propagator diagrams are shown in Fig. |6] where the 
doublon lines with squiggles represent the full doublon 
Green function to be obtained self-consistently and the 
thick single lines are the Fermion propagators. 

At this point, we make an additional approximation, 
and replace the k-dependent vertex functions Ak by mo- 
mentum averaged vertex functions y/ (A^) listed in Ta- 
ble [ij The basis of this approximation, is that within our 
resummation scheme, the vertex functions always occur 
in pairs with identical and largely arbitrary momentum 
indices, as can be seen from the self-consistent equation 
represented in Fig. [6] The self-consistent equation, there- 
fore, contains the product of this pair of vertex functions, 
and we replace this product by its momentum averaged 
value. 

Having replaced the momentum-dependent vertex 
functions by momentum-independent ones, we can re- 



place Green functions and self-energies by their momen- 
tum averaged counterparts. With this modification, the 
doublon self-energy is given by 



Erf(w) = zJ^C (w)-2zJ^ 



(13) 



T:^{lo)^zJ^C\lo)-2zJ^ r —S" {uj')g'^{Lo-u') 



— oo 

°° doj' 



gAio')S {uj-u:') (14) 



where S^, C, and S are the retarded doublon self-energy, 
Fermionic particle-particle and particle-hole propagators 
respectively [21 , and the primes ' and " denote the real 
and imaginary parts, respectively. The particle-particle 
and particle-hole propagators are given by 



5* {uj 
S\uj 

C"{uj 
C'{uj 



——Qf {^')Gf {uj'-u}) 





(15) 



J/ 

— gfiLu')[gf{Lu'^Lu) + gf{u;'+u;)] (16) 



— ^/f (w')^^f ('^ - w') 

TT J J 



(17) 



g"j{Lo')g){uj' + Lo)] (18) 



where gf{io) — X^k^/C*'^) momentum averaged 

Fermion Green function and the primes ' and " denote the 
real and imaginary parts. These equations, together with 
the equation for the doublon Green function, Eq. (12 1, 



define a system of self-consistent equations for the dou- 
blon self-energy. 

In this section we have made two approximations: (1) 
we replaced the momentum dependent vertex functions 
by momentum independent ones, and (2) we have left 
out a large number of diagrams with crossing Fermion 
lines (see Fig. [t] for some typical examples). To ver- 
ify these approximations, we have explicitly computed 
all diagrams up to 6*'* order in a Fermi Golden Rule 
calculation, which is free of these approximations (see 
Appendix [b] for details) . We find that the decay rate 
computed via Fermi Golden Rule matches very well with 
the resummation result. Further, within Fermi Golden 
Rule calculation we empirically verify that the contri- 
bution of crossed diagrams to the doublon self-energy is 
indeed negligible. Intuitively, the reason for this seems 
to be that the Fermion-doublon interaction vertex con- 
tains the factor 7p_q + 7k + 7q which changes sign as 
we sample momentum space. The non-crossing diagrams 
involve squares of this vertex function and do not change 
sign as we integrate over momentum coordinates. On 
the other hand, the crossing diagrams involve product 
of the vertices at different momenta and hence give a 
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FIG. 7: Typical diagrams for the doublon Green function 
that are not accounted for in the resummation procedure as 
they contain crossing Fermion lines. Here, double-lines stand 
for bare doublon propagators and thick single-lines the full 
(resummed) Fermion propagator. For reasons explained in 
the main text, these crossing diagrams do not contribute to 
the doublon decay, as the vertex factors are not paired and 
thus average to zero upon momentum integration. To see this, 
the momenta and a pair of vertices in the first, pretzel-like, 
diagram are labeled. Notice that the 7 vertex factors have 
different momentum labels, these would have been identical 
for the case of a non-crossed diagram. 




FIG. 8: Self-consistency equation for the Fermion propagator 
(top) and some typical diagrams that make up the full propa- 
gator (bottom). Thin lines indicate the bare propagator and 
thick lines the full (resummed) propagator. 



negligible contribution upon integrating over momentum 
coordinates. 



B. Fermion Self-Energy 

We now come back to the question of evaluating the 
Fermion Green function 



(19) 



where Ck = — >/7k — A* is the bare dispersion and S/(cl>) 
is the Fermion self-energy that arises due to interaction 
with other Fermions. To make progress, we begin by 
considering the non-crossing approximation. As before, 
for the case of the doublon self-energy, we are interested 
in the high frequency part of the Green functions, and 
therefore (in the non-crossing approximation) we are jus- 





FIG. 9: Typical crossed fermion diagrams that are missed 
by the resummation method. These types of diagrams are 
expected to strongly contribute to the Fermion self-energy at 
high frequencies and therefore to the doublon decay rate. 




FIG. 10; Typical type III diagrams that are missed by the 
resummation method. As explained in the main text, these 
diagrams are expected to give some contribution to the dou- 
blon self-energy, but their effect is not taken into account. 



tified in replacing the vertices by their momentum aver- 
aged counterparts as listed in first row of Table |Tj and 
then working with momentum averaged Green functions 
and self-energies. In the non-crossing approximation, the 
Fermion self-consistency equation is depicted diagram- 
matically in Fig. [8j where the thick Fermion lines repre- 
sent fully dressed Fermion Green functions that are be- 
ing determined self-consistently. The Fermion self-energy 
equations are given by 



S" {lo')G" [u - J) (20) 



+2zJl 



d^' 



g''{L0')S\LJ~u'). (21) 



Combining these self-energy equations with the definition 
of the Green function Eq. (19), we obtain a set of self- 



consistent equations for the Fermion Green function. 



C. Corrections Due to Diagrams Left Out 

In the resummation formalism we have missed three 
important classes of diagrams: type I diagrams, which 
correspond to doublon self-energy diagrams with cross- 
ing Fermion lines (examples depicted in Fig. [7| ; type II 
diagrams, which are Fermion self-energy diagrams with 
crossing Fermion lines (examples depicted in Fig. [9]); and 
type III diagrams, which are doublon self-energy dia- 
grams which are left out and are neither type I nor type 
II (examples depicted in Fig. 10). 

As mentioned earlier, we have empirically checked that 
type I diagrams do not contribute to the doublon self- 
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energy due to the lack of pairing of the Fermion-doublon 
vertex factors. However, there are no similar arguments 
for excluding type II or type III diagrams. We suppose 
that when a doublon emits a particle-hole pair, the parti- 
cle and hole are not coherent with each other, and there- 
fore, we make the approximation of dropping type III 
diagrams. However, each Fermion in the emitted pair 
still interacts with the Fermi sea, resulting in both non- 
crossing Fermion self-energy diagrams, that have already 
been taken care of, and type II diagrams which we shall 
try to estimate. 

Since we cannot evaluate all the type II diagrams ex- 
plicitly, we proceed to approximate their effect on the 
Fermion self-energy in the following way: 

(a) We assume that at a given frequency w, the lead- 
ing contribution to the imaginary part of self-energy 
Im Y.f{uj) comes from diagrams of a definite order no(w), 
as diagrams of lower order do not have enough particle- 
hole pairs to absorb w and diagrams of higher order 
are suppressed by additional powers of J/uj. We ex- 
pect no(a;) to scale linearly with uj as the main contri- 
bution to the spectral function at uj comes from exciting 
~ uj/eo particle-hole pairs, where eq is the typical energy 
of particle-hole pairs. 

(b) To determine n(){uj), we keep the Fermion- Fermion 
vertex energy scale Ji as a free parameter, and calculate 
no{uj) from the logarithmic derivative 



no{uj) 



1 dlogI]/(a;) 



2 dlogJi 



(22) 



This relation is exact if only one order of diagrams con- 
tribute at given energy; for the case of different orders 
contributing to self-energy, this gives a number close to 
the order with leading contribution, rio(w), obtained 



from the resummed self-energy, is plotted in Fig. (111. 
The best fit for this graph is no{uj) = a;/(5.85J) — 1/27 

(c) We then compute the ratio of the total number of 
possible n}^ order Fermion self-energy diagrams to the 
number of n^^ order diagrams included in the resumma- 
tion scheme, R{n). R{n) can then be interpolated to form 
a function of the continuous variable n. See Appendix [C| 
for details of computing this ratio. 

(d) In the final step, we rescale the imaginary part of 
the Fermion self-energy by R{nf){uj)) to obtain a better 
approximation including effects of missed diagrams 



(23) 



Here, we are making an assumption: the amplitude of 
the Fermion self-energy diagram only depends on its or- 
der in perturbation theory and not on the details of the 
structure of the diagram. Modulo the contribution of 
the type HI diagrams, this approximation should over- 
estimate the decay rate as the crossing diagrams usually 
contribute less than the non-crossing diagrams due to the 
momentum sums involved. 

To complete the calculation of the doublon self-energy, 
we use the Fermion Green function to construct the 
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FIG. 11: Order with largest contribution to the Fermion self- 
energy no{u}) as a function of the frequency lu. The solid line 
represents the best linear fit for the high frequency data. 



particle-particle and particle-hole propagators Eqs. (15 



18 1, which appear in the self-energy equations for the 
doublon Eqs. (pinil. 



V. THEORETICAL RESULTS AND 
COMPARISON WITH EXPERIMENTS 

In this section we look at the theoretical results of the 
doublon lifetime calculation and compare them with ex- 
perimental results. We start by summarizing the method 
of calculation, which will help in establishing different ap- 
proximation schemes. We then discuss the results from 
different schemes and their comparison with experiments. 

The calculation of the decay rate via the resumma- 
tion technique has two important steps. The first one 
is the evaluation of the Fermion Green's functions which 
are used to compute the particle-particle and particle- 
hole propagators. The second one is the evaluation of 
the doublon self-energy, which uses these propagators. 
As mentioned before, a non-crossing approximation for 
the doublon self-energy yields good results. The crossing 
diagrams give negligible contribution as the vertex func- 
tions which oscillate with momenta kills the momentum 
averages. We also note that there is a set of doublon 
self-energy diagrams (the type HI diagrams) which we 
neglect in our calculation. 

Our approximations are then related to different 
ways of evaluating the Fermion propagators. We con- 
sider three different approximations: (i) Non-interacting 
Fermions; in this case we use the free Fermion propaga- 
tors with a band dispersion. One way of looking at this 
approximation is to set Ji = 0. (ii) Non-crossing ap- 
proximation for interacting Fermions; in this case we set 
Ji = J but use only non-crossing diagrams to evaluate 
the Fermion propagators, (iii) Modified self-energy for in- 
teracting Fermions; in this case we modify the self-energy 
of the interacting Fermions obtained by non-crossing ap- 
proximation to take into account Fermion self-energy di- 
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FIG. 12: Fermion spectral functions in different approxima- 
tions: the free Fermion spectral function {^AS^^\liS)']\ the pro- 
jected Fermion spectral function obtained as the result of the 
resummation procedure (A(a;)); the projected Fermion spec- 
tral function including corrections for missing diagrams in the 
resummation procedure (corrected A(lS)'). The linear slope 
at high energies on a semi-logarithmic scale shows the expo- 
nential transfer of spectral weight due to projection induced 
interactions. 



agrams missed in the resummation. The modification 
procedure is detailed in the previous Section. 

We plot the spectral function of the Fermions, A{bj) 



— (l/7r)lmtjy (w), for the three approximations in Fig. 12 
In the non-interacting case, this is simply the density o 
states in a cubic lattice and the spectral weight is zero 
outside the band. In the non-crossing approximation, we 
see that there is a transfer of spectral weight from low 
energies to an exponential tail at high energies, which 
reflects the fact that interaction induced by projection 
leads to the possibility of creating a high energy Fermion, 
which can reduce its energy by creating particle-hole 
pairs. This is an important qualitative change that af- 
fects the physics of doublon decay in a fundamental way. 
The interacting Fermion approximation allows two dis- 
tinct decay processes : (a) creation of several low energy 
(w ~ 22; J) particle-hole pairs and (b) creation of a high 
energy particle-hole pair which then decays into a shower 
of low energy particle- hole pairs. The second process is 
forbidden for non-interacting Fermions. Finally, in the 
modified self-energy approximation, we include more pro- 
cesses to create particle-hole pairs and hence there is a 
larger shift of spectral weight to higher energies, as ev- 
idenced by the slower decay of the tail. This enhances 
the importance of the (b) channel for decay. 

In the second step we use the Fermion propagator ob- 
tained in step one to self-consistently compute the dou- 
blon self-energy. The imaginary part of the doublon self- 
energy for various U /QJ ratios is depicted in Fig. 13 The 
main features are a pair of peaks, one occurring at small 
frequencies, and another at high frequencies. As there 
are no excitations in the Fermi system in the initial state, 
for frequencies uj < U a, nonzero value of ImE(;(a;) cor- 
responds directly to the rate of doublon decay. At low 



FIG. 13: Doublon self-energy (in the modified Fermion self- 
energy approximation) for various values of U/6J. 
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FIG. 14: Doublon decay time as a function of U /QJ . The 
blue circles are the experimental data (cf. Fig. |4|. The lines 
represent theoretical results from resummation with different 
levels of sophistication from non-interacting Fermions (red 
dashed line) to the non-crossing approximation with inter- 
acting Fermions (green dot-dashed line) to the modified self- 
energy approximation (purple solid line). 



frequencies, the doublon is far from its mass shell and 
rapidly decays into a pair of particles. As the frequency 
increases more and more particle-hole pairs are required 
to absorb the doublon energy resulting in the exponential 
decrease in ImSrf(w). As uj surpasses U, a new contri- 
bution to the imaginary part of the doublon self-energy 
arises from processes where the doublon can scatter into 
a lower energy state closer to the mass shell by releasing 
the excess energy in the form of a few particle-hole ex- 
citations. This scattering process is responsible for the 
high frequency peak in lmS,d{^), that starts growing at 
u! = U. As we are interested in the decay of a doublon 
on the mass shell, we read it from ImSdiU), which cor- 
responds to the smallest value of lml](i(aj) between the 
two peaks. 

In Fig. [T4j we plot the experimentally obtained de- 
cay time together with the theoretical estimates from 
the three different approximations mentioned earlier. We 
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proceed in the order of sophistication, starting from the 
non-interacting Fermion case. We see that the decay 
time obtained with non-interacting Fermions (Ji — 0) 
via resummation of doublon self-energy diagrams is much 
longer than the experimentally obtained one. Setting 
Ji = J, and using non-crossing diagrams for Fermion 
self-energy, we obtain a decay time that is a closer match 
to the experimental data, but is still too long. Next, we 
take care of the corrections to the imaginary part of the 
Fermion self-energy from crossing diagrams and find a 
reasonable match with experiments. 

Finally, we want to comment on the remaining free 
parameters in our calculation. The chemical potential of 
the Fermions, which determine the hole density, is a free 
parameter, which can in principle be determined from an 
equilibrium theory of a strongly interacting doped Hub- 
bard model. Since there is no consensus about the theory 
of the doped Hubbard model, we prefer to keep it as a 
free parameter. We vary it within the plausible range 
of 0.25 J to (— 0.3J) to see how sensitive our results are 
to the choice of this parameter. The dispersion in the 
lifetime is then plotted as the shaded region in Fig. [4] 
We see that we find good quantitative agreement with 
the experiments in the slope of the lifetime curve, i.e. for 
the co-efficient a in the exponent of the scaling function. 
The agreement in the prefactor C is also fair, but this 
quantity is sensitive to the choice of the free parameter 
in our calculation. 



VI. CONCLUDING REMARKS 

We have studied the decay of artificially created dou- 
ble occupancies in the repulsive Fermi-Hubbard model in 
the presence of a background compressible state. The sit- 
uation is experimentally realized by creating double oc- 
cupancies and corresponding holes on top of a half-filled 
system via optical lattice modulation. Experimentally it 
is found that the decay time of the doublons scales ex- 
ponentially with U/J. We can understand the observed 
scaling in terms of the fact that in order to decay the 
doublon has to distribute its energy (~ U) among ^ U / J 
particle-hole excitations. We have developed a detailed 
theoretical description of this process using diagrammatic 
resummation techniques. Although the scaling form can 
be understood from a simple energy conservation argu- 
ment, we find that the co-efficient in the exponent de- 
pends substantially on the strong interaction between the 
background Fermions. After taking into account the ef- 
fects of these strong interactions, we find quantitatively 
fair agreement between theory and experimental results. 

The exponentially large lifetime of the doublons has 
serious implications for use of cold atom systems to sim- 
ulate the equilibrium properties of the Hubbard model 
at large values of U / J . Typically, in cold atom exper- 
iments, the strong interaction regime of the Hubbard 
model is accessed by cooling the atoms in a weakly in- 
teracting state and then tuning either the optical poten- 



tial or the magnetic field to change U/J. The lifetime 
of the doublons constrains the maximum sweep rate of 
these Hamiltonian parameters under which thermal equi- 
librium is maintained. As one goes towards larger C//J, 
the sweep rates need to be exponentially slow to maintain 
thermodynamic adiabaticity. Given intrinsic constraints 
like lifetime of a sample, this would restrict the values 
of U I J for which the simulation of Hubbard model in 
thermal equilibrium can be achieved. 

However, this also opens up the possibility of studying 
non-equilibrium dynamics of the Hubbard model, which 
may contain interesting and new physics. In addition, 
the long lifetime of the doublons also leads to the possi- 
bility of observing metastable states with finite density of 
doublons. An intriguing scenario is observing r\ pairing 
of doublons and holes [22j . 

Finally, we point out that similar phenomena may be 
relevant to the issues of equilibration in the Bosonic Hub- 
bard model. In a recent paper [53] C. Chin's group 
observed the equilibration of the density distribution of 
Bosonic atoms in a two dimensional optical lattice after 
the lattice potential was ramped up. As the system re- 
laxed toward equilibrium, the center of the trap heated 
up, which required the increase in the number of dou- 
blons. The slow relaxation timescale observed in experi- 
ments may be a reflection of the "dual" problem to the 
one we discussed in this paper: slow rate of formation of 
doublons from a a state containing only singly occupied 
sites and holes. 
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Appendix A: Model 

In this appendix we derive the model we use to describe 
doublon decay in the background of a projected Fermi 
sea. We begin with the Fermi-Hubbard model 

i^FH = -'^Y. + ^ II (Al) 

(ij)cr i 

where the first term describes the hopping of fermions 
and the second term the on-site repulsive interaction. We 
are interested in the case U ^ J, where we expect dou- 
blons to be meta-stable particles. Therefore, our goal is 
to decouple the doublon sector from the sector of singles. 
We do this by projecting out double occupancies from 
the singles sector, and introducing doublon creation and 
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annihilation operators d] and di to take their place. We, 
proceed in two steps, first we use projection operators to 
separate the terms in the Fermi-Hubbard Hamiltonian 
that preserve the number of doublons from those that 
change it: 

Hpii = Ho + H+i + H^i, (A2) 
where -ffo preserves the number of doublons 

i 

and i?±i increases/decreases it by one 

H-i = 51 (1 - n^MaCjainjs), (A5) 

where riia = cl^c^„ and a indicates spin opposite to cr. 
In the second step, we replace double occupancies by the 
corresponding doublon operators. Thus we have 

-JY, 4d^c,,cl + uY,nt (A6) 



and 



(A7) 



H+i = -J Y <7dl{CjgCia){l ~ n,g) 
H-1 ^-JT. - n,<^){clc]g)d,, (A8) 

(u>o" 

where nf = d^d^. Thus far, we have obtained an expres- 
sion for the Fermi-Hubbard Hamiltonian that incorpo- 
rates doublon operators. This Hamiltonian was specifi- 
cally derived in such a way as to avoid creation of spu- 
rious states (e.g. a doublon and a single fermion on the 
same site) by the use projection operators. As a result, 
we do not need to supplement it with a constraint equa- 
tion. 

Now we can separate the terms in the Hamiltonian 
based on which sectors they connect. The Fermion- 
Fermion term arises from terms in Hq that connect the 
projected sector and is given by 



Likewise, the Doublon repulsion term also arises from Hq 
and is given by 



Ha^uY^nf. 



(AlO) 



The remaining terms connect the Fermion and Doublon 
sectors and are 



Hjd = H+i +H^i + jY 



(All) 



(1 - n,g)n'^ + nf{l - n,g) + d]d] 4c^.,, (A12) 



where we have dropped the term that is nonzero in the 
presence of a pair of doublons as we are assuming that 
there is at most one doublon. To complete the model, 
we drop terms that result in Feynman vertices with more 
than two incoming and two outgoing propagators. We 
have verified, numerically, that these diagrams do not 
significantly contribute to the doublon decay rate. 



Appendix B: Checks on Approximations through 
Fermi Golden Rule Calculation 



In this appendix, we compute the doublon decay rate 
for the case of non-interacting Fermions (i.e., we disre- 
gard Hp part of the Hamiltonian (Isl)). We treat H^ = 



Hj + Hd as the base Hamiltonian, and Hfd as the per- 
turbation Hamiltonian, and evaluate the decay rate, via 
the Golden Rule, to very high order in Hfd using Monte 
Carlo integration. The objective of this appendix is to 
test the approximations made in the resummation tech- 
nique of Section |IV| on a simplified Hamiltonian. In par- 
ticular, we empirically verify that (1) we may ignore the 
crossing diagrams in doublon self-energy and (2) we can 
use momentum averaged Green functions to compute the 
decay rates. We begin by laying out the formalism, and 
then list the results of Monte Carlo integration of decay 
rates. 



1. Formalism 

Our goal is to compute the transition rate from the 
starting configuration composed of a single doublon in 
a Fermi sea at finite temperature to the final configura- 
tion composed of the initial Fermi sea with the doublon 
converted into a pair of single particles and a number of 
particle-hole excitations. The Fermi Golden rule states 
that the decay rate is given by 



27r 



(Bl) 



14 



where the matrix clement can be expressed in ordinary 
perturbation theory via 



)(Sn-l\Hfd\Sn-2)--{si\Hfd\i) 



{E^ — Es^){Ei — Es2)...{Ei — -Es„_J 



(B2) 

Here, the sum goes over aU intermediate states Si, with 
energy Eg. , and n is the order of perturbation theory. 
In this perturbation theory, the action of Hfd (except 
for the final matrix element {f\Hfd\sn-i)) is to create 
particlc-holc pairs. In principle, we may be able to con- 
nect the initial state to the final state via other processes, 
e.g. doublon^particle-particle— >-doublon, however, these 
process lead to decay at higher order in perturbation the- 
ory, and thus we ignore them. 

We label the initial state by the momentum of the dou- 
blon p: 



|i) = |-;p)=4|FS). 



(B3) 



Likewise, we label the final state via a set of momenta 
for the up (down) spin particles ki^(^) and the up (down) 
spin holes g^j^j.): 



\f) = 



(B4) 



X «,t^«n,t) - |FS), (B5) 



where n-^(i) counts the number of spin up (down) 
particle-hole pairs created (n-|- + ni + l = n). 

The intermediate states arc composed of a donblon 
and 1,2, 3,..., n — 1 fermion-hole pairs. Using Hf^, we 
can write the matrix element as 



sig(perm) 



T\-,P) 

{f\Hi\{ai,ki,qi), (a„_i, fc„_i, g„_i);p„ 



i)...{{ai,ki,qiy,pi\Hi\p) 



permutations 



i^ki + ••• + ^kn-i ^91 



where fcj, qj, py stand for the particle, hole, and dou- 
blon momenta, respectively, and di indicates the spin 

of the i-th particlc-holc pair. The sum runs over all 
intermediate states that lead to the final state |/). 
That is, we must sum over all permutations of as- 



signed values to [iji, ki, qj) from the list {k 



{fci,^,...,fc„^+i}, {qi,^, ...,(7„J, {qi^i,...,qnj. Within 
this labeling scheme, the doublon momenta in the inter- 
mediate states py, and the hole momentum in the final 
state, are chosen automatically by momentum conserva- 
tion. We take care of the Fcrmionic anti-commutation 
relations with sig(perm), which stands for the signature 
of the permutation, and is ±1 for even/odd permutations 



(B6) 
(B7) 



r 



of momenta. 

To obtain the decay rate, we trace over the final states, 
order by order in perturbation theory, 



r(p) = ^r"(p). 



(B8) 



n=0 



At each order we trace over the number of up- and down- 
spin particlc-holc pairs, and the corresponding momenta 
of particles and holes that make up the final state. The 
decay rate at n-th order is the given by the expression 



r"(p) 



27r 



[5fci,i-...5fc„^+i,-|-] [5fci4...5fc„^+i4] [5gi,i-...%nt,t] [%,4.-%n44] 



nf+ni + l=n 



(nt + + l)!(nt)!(n4)! 

5iU-Ef)5(p-Y,k + T.'i) 



A:i_^...fc„^+i,-|-, ki^^...kn_^_+i,i, 



T\-,P) 



(B9) 



where dk stands for f {k) (P k / (2^)^ , dq for (1 — /(g)) (i^g/(27r)^, and /(fc) is the Fermi function. The de- 
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FIG. 15: Decay rate as a function of the order of the per- 
turbation theory computed using Fermi Golden rule. Largest 
decay rate corresponds to most important order. 
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FIG. 16: Comparison of the resummation method and various 
Golden Rule approximations for calculating the dependence 
of the Doublon decay time on the interaction strength (7/6 J 
(with non-interacting Fermions). 



nominator in the integral takes care of the fact that inter- 
changing a pair of momentum labels does not change the 
final state, = ^(fci,^) ... ^(fcnt+i^T) + ^(^i,;) + - + 
C(A:n^-M,;) -^(9i,t) - - -S(9nt,T) -^('Zi.i) - - -C(9n^4)) 
is the final state energy, and the second (5 function takes 
care of momentum conservation. 

We explicitly evaluate the 3^" dimensional integral in 
Eq. (B9| numerically via Monte Carlo integration. To 



perform this integration, we replace the 6 function of en- 
ergy, which defines a hypersurface in momentum space - 
a volume of of measure zero, by the top hat function. We 
also use important sampling to speed up integration by 
biasing our selection so that we pick particle-hole pairs 
with holes in the Fermi sea and particles outside of it. 
The main numerical constraint on the speed of integra- 
tion comes from evaluating the (n^ -t- \)\(n\^ 4- l)!ri^!n^! 
permutations over the intermediate states, which be- 
comes rather expansive for n > 6. 



Results 



We begin by verifying that the perturbation theory 
in Hfd does indeed converge. That is, for fixed J7/6J, 
does r("){p) decrease sufficiently fast as n increases? We 
know that for n < U/12J, r'^"^(p) = 0, as not enough 
particle-hole pairs are formed to carry away the energy 
of a doublon. When n ~ U/12J, in order to satisfy 
energy conservation, particles created in the decay must 
have momentum in vicinity of the band maximum near 
(tt, tt, tt) and holes in the vicinity of the band minimum 
at (0,0,0). Therefore, for n ^ U/12J the volume of 
the momentum space being integrated is very small, but 
this volume increases quickly as n grows. As a result, 
we expect that the F^"^ (p) will increase with n for small 
n. On the other hand, at high orders the decay rate is 
suppressed by a high powers of the small parameter J/U. 
Thus, we expect F(")(p) to have a maximum for some 
intermediate value of n close to, but somewhat larger 
than U/12J. 



In Fig. 



15 



we plot F'^")(p) as a function of n for various 
values of U/6J. In all cases, computations have been 
performed at T = and fi — (corresponding to one 
particle per two sites). As expected, in all cases, we see 
a clear peak in F(")(p) at n U/12J + 2. 

Having verified the convergence of the high order per- 
turbation expansion, we move on to empirically verify 
whether we can ignore crossing diagrams, at least for 
the case of free Fermions. In order to perform this com- 
parison we compute the total decay rate as a function of 
U/ 6t using both Monte Carlo integration of Eq. ( B9 ) (in- 
corporates all possible diagrams), as well as the resum- 
mation of the non-crossing diagrams given by Eq. ( 14 ) 
with bare Fermion Green functions used to compute C{uj) 
and S{lu). We perform two additional tests using Monte 
Carlo integration: (1) We calculate the decay rate with 
Bosonic instead of Fermionic signs for closed Fermion 
loops; (2) We keep only the diagonal terms, i.e. we re- 



place 



^ perm ■ 



^ perm I 



which corresponds to the 



order-by-order summation of non-crossing diagrams, but 
without momentum averaging of the resummation ap- 
proach. The results of these four types of calculations 
are plotted in Fig.[T6j for T = and /i = 0. There is very 
good agreement between all four cases, confirming that 
crossing diagrams may indeed be dropped as explained 
in subsection IIVAI 



n=2 




FIG. 17: All distinct tree diagrams with one vertex (left) and 
two vertices (right). 
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FIG. 18: Dependence of the number of distinct tree diagrams 
on the number of nodes in the tree. 



Appendix C: Diagram Counting 

In this appendix we describe the procedure for count- 
ing the total number of distinct, spin-labeled Fermion 
self-energy diagrams at a given order Qediin) and the 
number of non-crossed spin-labeled Fermion self-energy 
diagrams Qnc{fT')- We remind the reader that (3aii(?^) and 
Qnc('^) correspond to diagrams with 2n vertices. For high 
uj, Yij:{uj) is dominated by diagrams with maximal num- 
ber of particle and hole lines in the middle, as these max- 
imize the energy that is being transferred to the particle- 
hole pairs being created. In fact, the range in uj over 
which Sy(i^) is nonzero is proportional to the number 
of particle- and hole-lines in the middle of the diagram. 
Therefore, to simplify the counting, we only count dia- 
grams that have the maximal number (2n+l) of particle- 
and hole-lines going across the middle of the diagram. 

To count the number of diagrams at given n, we first 
construct all distinct tree diagrams (without spin labels) 
that have a single particle going in, n -|- 1 particles and 
n holes going out and n vertices of the type given in 
first row of Table |l] In Fig. [TTj we show all such tree 



diagrams for n = 1 and n = 2. In Fig. 18 we show 
how the number of distinct trees scales with n. Next, we 
construct the set of all the possible self-energy diagrams 
by taking a pair of tree diagrams, reversing all the arrows 
in one of them, and gluing them together. When we 
count the total number of diagrams, we glue together 
particle-particle lines and hole-hole lines in all pairs of 
trees at the given order, in all possible ways. On the other 
hand, when counting the number of diagrams produced 
by the non-crossing approximation, we only glue together 
trees with their mirror image. Finally, we spin label the 
resulting diagrams, and remove all duplicate diagrams, 
to obtain Qa,ii{n) and Qnc{n)- 

We assume that the ratio (3au("-)/Qnc("-) scales like 
^ e"". We use this assumption to extrapolate the ratio 
for non-integer values of n and for large values of n > 
4. We plot the ratio of Qaii{'n)/Qnc{'n), along with the 
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FIG. 19: Correction ratio as a function of the order of the 
diagram. 



extrapolated curve that we use in rescaling the Fermion 
self-energy, in Fig. [19} 
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